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O ■ ABSTRACT 

o ■ 

(N 

^ I We devise a method to measure the abundance of satellite halos in gravitational 

lens galaxies, and apply our method to a sample of 7 lens systems. After using 
Monte Carlo simulations to verify the method, we find that substructure comprises 
fsat = 0.02 (median, 0.006 < fsat < 0.07 at 90% confidence) of the mass of typical lens 
CS| ■ galaxies, in excellent agreement with predictions of CDM simulations. We estimate a 

characteristic critical radius for the satellites of O'.'OOOl < b < 0'.'006 (90% confidence), 
i For a dn/dM oc M~^'^ (Miow < M < Mhigh) satellite mass function, the critical radius 

' provides an estimate that the upper mass limit is 10^ Mq <^ Mhigh ^ 10^ M©. Our 

measurement confirms a generic prediction of CDM models, and may obviate the need 
to invoke alternatives to CDM like warm dark matter or self- interacting dark matter. 

^: 

I ' Subject headings: cosmology: theory - galaxies: formation - gravitational lensing 

■ large-scale structure of universe - dark matter 

•I— I . 1. Introduction 

X 

■ A discrepancy between the number of satellite halos expected from CDM simulations and 
the observed numbers of Galactic satellite galaxies is part of the prosecution's case for a crisis in 
the CDM scenario for structure formation (e.g. Kauffmann 1993, Moore et al. 1999| , Klypin et 



al. 1999| ). Suggested solutions range from the mundane, such as the inhibition of star formation 



in the satellites by photoionization (e.g. Klypin et al. 1999, Bullock, Kravtsov & Weinberg 200C| ), 



to the exotic, such as the disruption of the satellites by self-interacting dark matter (e.g. Spergel 
& Steinhardt 2000| ) or changes in the power spectrum (e.g. Bode, Ostriker & Turok pOOl , Colin, 



Avila- Reese & Valenzuela 2000 ). The satellite crisis must also be closely related to the more 
general problem that the number of low luminosity galaxies diverges only as 1/L ~ 1 /M while 
the number of CDM halos diverges as ~ 1/M^, implying that the probability of forming a visible 



galaxy in a low mass halo must diminish as ~ M (e.g. Scoccimarro et al. [2001, Kochanek 2001 



Chiu, Gnedin & Ostriker 2001| ). In principle, the measured abundances of satellite halos should 



provide a strong test of the CDM scenario, but because the satellites used as evidence for a 
problem have low luminosities and (in many cases) low surface brightness, it is difficult to apply 
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this test to any galaxy besides our own. Moreover, the test only considers the numbers of satellites 
with detectable optical emission, which is at best a lower bound on the number of CDM halos. 

Gravitational lensing is the only probe which avoids both of these limitations, as was already 
noted by Moore et al. ( 1999 ). First, the test can be applied to many lens systems spanning a 
range of redshifts and physical properties. Second, because lensing phenomena couple directly 
to mass, lenses are sensitive to both luminous and dark substructures in CDM halos. Mao & 



Schneider ( 1998 ) pointed out that the anomalous image flux ratios observed in several lenses, 
particularly B1422+231, could be explained by substructures such as low mass satellites in the 
primary lens galaxy. The primary lens magnifies the perturbations from the substructure, making 
the brightest images particularly susceptible to the effects of substructure. Recently, Metcalf 
& Madau ( |200lD quantified the effects of CDM satellit es using simulations and found that the 



effects should be readily detected, and Chiba ( 2001 ) demonstrated that plausible CDM satellite 
distributions could explain the anomalous flux ratios in B1422+231 and PG1115+080. Detailed 
studies of B1422+231 (Keeton |2002 , Bradac et al. 2002| ) find that the observed perturbations 
require substructure with mass scales comparable to CDM substructure (^ 10® Mq) rather than 
stellar microlensing, and Metcalf &: Zhao ( p002| ) have shown that the anomolous flux ratios cannot 
be reproduced in a large family of smooth potentials for the primary lens. 

The missing link is an approach for analyzing the gravitational lens data to estimate the 
properties of the satellite population. In this paper we develop such an analysis method and apply 
it to a sample of 7 lenses to estimate the surface density and characteristic mass of the perturbing 
satellites. We focused on analyzing four-image radio lenses because using the radio lenses 
eliminates the problem of dust extinction, and minimizes the problems from stellar microlensing 



due to the relatively large source size (see Koopmans & de Bruyn 2000). We analyzed the lenses 
MG0414+0534 (Hewitt et al. |199^ ), B0712+472 (Jackson et al. |l998D , PG1115+080 (Weymann 
et al. |198C| ), B1422+231 (Patnaik et al. |99|), B1608+656 (Fassnacht et al. |l99^ ), B1933+503 
(Sykes et al. |199^ ) and B2045+265 (Fassnacht et al. |l999| ). Of these 7 four-image lenses, 6 
show anomalous flux ratios which might be due to the effects of substructure. We develop our 
formalism, characterize our model for the satellite distribution and test our analysis methods in 
§§. We apply it to the lens sample in §|3[ In §Q we review our conclusions and their limitations and 
then outline the observations needed to improve them. 



2. Analyzing the Effects of Substructure on Gravitational Lenses 

In this section we outline our mathematical approach to analyzing the lenses to determine the 
properties of the substructure (§2.1) and the physical model we use for the satellites composing 
the substructure (§2.2). In §2.3 we discuss the relationship between our model and the physical 
properties of the substructure such as its fractional surface density, mass and velocity scales and 
linear sizes. In §2.4 we outline our Monte Carlo models and test the analysis method. 
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2.1. A Linearized Approach to Analyzing Substructure 

Unlike the primary lens galaxy, which we can observe directly to determine its position and 
optical properties (e.g. Lehar et al. 200C , Kochanek et al. 2000| ), we can detect substructure only 



through its effects on the positions and fluxes of the lensed images. This means that estimates for 
the properties of the substructure will be difficult to separate from the properties of the primary 
lens (the "macro" model to adopt the language of the quasar microlensing literature) because 
many perturbations will be degenerate with changes in the macro model. For this reason, Mao & 



Schneider (1998) focused on merging image pairs where macro models must generically predict 
similar image fluxes but the observations sometimes find very different fluxes. We will instead 
allow the macro model to compensate for the effects of substructure as part of our analysis. If 
we confine our analysis to typical four-image (two-image) lenses we have 14 (8) constraints for 
determining the 10 parameters of a realistic macro lens model. For a four-image lens, the macro 
model is overconstrained and we can attempt to estimate the properties of the substructure. 
Because we have typically found that it is relatively easy to fit image positions and hard to fit flux 
ratios using standard lens models, we expect the deflection perturbations from substructure to be 
small or degenerate with the parameters of the macro model, and the magnification perturbations 
to be larger and non-degenerate. 

We model the lens by combining a macro lensing potential </>(x, p) defined by a set of 
parameters p, with a localized perturbing potential for each image (50j(x). For later notational 
simplicity, the source position and flux are considered part of the parameter vector p. The time 
delay surface near image i is 

T = ^(U - X)2 - 0(x,p) - 50i(x) = To(x,p) - 50i(x), (1) 

and we find images at solutions of Vr = with an inverse magnification tensor of M^^ = VVr 
and magnification M = |M| (see Schneider, Ehlers & Falco 1992| ). 



We assume the substructure produces small perturbations, so we can simplify its effects 
by expanding the lens equations as a linear perturbations to a macro model for each image i. 
We would like to linearize the equations so that the process of adjusting the macro models to 
compensate for the effects of substructure can be done rapidly. For macro parameters po we find 
images at positions ^f^^ with magnification tensors m|'^^ at the solutions to Vto(x.''\po) = 0. 
Expanding the lens equations about these solutions, the perturbed image positions are 

xf ) = xf ) + Mf ) • (<5x, - Ap • Ci) (2) 

where 5xj = V5(j)i is the deflection produced by the substructure and 

^ dVro 



dp 



(3) 



evaluated at x = x.f'^ and p = po is the change in the macro model deflections produced by a 
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where 



small change Ap = p — po in the macro model parameters. The perturbed image magnification is 

= M^°> + — — . + — • Ap + — • Ax (4) 

dSM dp dx ^ ^ 

which becomes 

jVf (1) = m(°) (1 + Srui + Ap • + S-Xi ■ Ej) (5) 

Srrii = -tr (mSM^^^ for ^M"^ = -VV(5(^i, (6) 

is the perturbation in the magnification due to the effects of the substructure on the magnification 
tensor, 

^, 1 fdM dM \ 

D = — — • M • C (7) 

is the perturbation due to changing the lens parameters, and 

E = — — -M (8) 

M dx ^ ^ 

is the perturbation due to the effects of substructure on the deflections. The functions 5mi, D- and 
E are all evaluated at x = xf^^ and p = Po- The flux of image i is /j = /sMj = + A/)Mj, 

which we can linearize as 

fi'^ = (1 + Srui + Ap • + ■ Ei) (9) 

where Ap • Dj = Ap • + A/ and the fractional change in the source flux A/ is considered one 
of the model parameters p. 

We detect substructure as residuals in fits to the lens parameters which cannot be modeled 
by the macro potential. If we use a statistic, the fit statistic for the image positions is 



1=1 



where (To ~ C'COS is the uncertainty in the observed image positions x"***. The fit statistic for the 
image fiuxes is 

A-/ - E f ~ ^''^^^ + (5m, + Ap • D + (5x, • E,) y ^^^^ 

i=l \ ^f'"'' ) 

for observed fluxes jf'^ and flux uncertainties a f^i ~ The lens position, if measured, 

is constrained by = (x"^'' — y^f^ — IS.pi)'^ /af where Ap; represents the perturbations to the 
lens position and ai ~ (/.'OOS is the uncertainty in the observed position x°''^. Because all the 
terms entering the fit statistic depend only linearly on the Ap, the fit statistic is a quadratic 
of the form = Xo + 2Ap • I + Ap • J • Ap which is minimized for Ap = — • I with value 
Xmin ^Xo^^-'J"^'-"-- If macro model were held fixed, the goodness of fit would be Xo) 
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combining the effects of the substructure with the differences between the observations and the 
initial macro model. After adjusting the macro model, the correction A^^ = I • J^^ • I represents 
the ability of the macro model to mimic the effects of substructure. 

In order to test the CDM predictions for substructure we need to estimate the statistical 
properties of the substructure rather than discuss the evidence for perturbations from substructure 
in individual lenses. For a statistical model of the substructure characterized by parameters s, the 
probability of the ath substructure realization for lens j, 6aj, is P{6aj\s). Given a concrete model 
for the substructure 6aj we can compute the likelihood for fitting the data as the likelihood of 
obtaining the statistic we find after reoptimizing the macro lens model given the substructure 
realization, P{Dj\6aj)- The Bayesian probability of the model parameters given the data for 
j = 1 ■ ■ ■ N lenses is 

P(s, 6o,j\Dj) a P{s)Uf^^P{Dj\6^j)P{6aj\s) (12) 

where P{s) is the prior probability distribution of the s. As in all Bayesian probabilities, the 
expression summed over all variables is normalized to unity. In estimating the statistical properties 
of the substructure, we are not interested in the likelihoods of the individual realizations, but in 
the marginalized distribution 

M 

Pis\Dj) oc P(s)nf=i PiDj\Sc.j)Pi5aj\s) (13) 

a=l 

where we sum over the a = 1 • • • M Monte Carlo realizations of the substructure for each lens. 
Typically we used M = 10^ realizations. As we vary the statistical properties of the substructure 
s, the fraction of the realizations 5aj which significantly improve the goodness of fit varies. These 
changes in the fraction of realizations which improve the fit relative to the macro model alone 
allow us to estimate the parameters s describing the substructure. 

In our final analysis we used random realizations of perturbing satellites to estimate the 
perturbations. It is worth mentioning, however, that the problem can be fully linearized if we 
use a Gaussian model for the perturbations. If k = {Smi,5xi} is a d-dimensional vector of the 
perturbation variables, and they have a covariance matrix = (k^k), then the Gaussian model 
for the probability distribution of the perturbations is 

P(k) = |V|i/2(2vr)-'^/2 exp (^-^k • V • k^ . (14) 

The matrix V is proportional to the inverse of the satellite surface density, so by combining 
eqns. (|To|), ( |Tll) ( [l^ and (IT^), the marginalizing integrals over the shifts in the lens model and 
the distribution of satellite realizations can be done analytically using standard methods for linear 
algebra and Gaussian integrals to leave an expression depending only on the statistical properties 
of the substructure. We did not use this for our actual analysis because it was relatively easy 
to perform the necessary Monte Carlo realizations and integrals needed to reproduce the true 
probability distribution for the substructure and its correlations, but we did use it as an internal 
check on our results. We mention it here because other studies may find it to be of similar utility. 
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2.2. Substructure Models 

Before applying the above formahsm to observed lens systems, it is necessary to calculate 
the expected amplitude of perturbations to lensed images from CDM substructure. This involves 
calculating the expected levels of astrometric shifts in the image positions, and the RMS 
fluctuations in the local convergence and shear. These quantities are then magnified by the local 
magnification tensor of the smooth macro-model. 

We model CDM substructure by randomly laying down subclumps of surface density. 
The substructures seen in CDM simulations appear to have mass profiles consistent with the 
'universal' NFW profile, however for simplicity we will treat them as pseudo-Jaffe models (density 
p oc r~^(r-^ + a?)~^, see Munoz, Kochanek & Keeton 2001), with convergence, the surface density 
in units of the critical surface density, 



S b 
Kir) = — = — 



1 



r (r^ + a^)-*-/^ 



(15) 



where the critical surface density for lensing is Sc = c^Dqs / ^t^GDolDls- Here, 6 is a length scale 
similar to the Einstein radius of the subclump, and o is a tidal or break radius. Note that 6, r and 
o are angular lengths, which are related to physical sizes by multiplication by the distance to the 
lens Dql- The total mass of a clump is M = vrftaSc where Sc is the critical density in angular 
units. If the surface mass density of the perturbers is S, the number of perturbers per unit area is 
N = (S/Sc)/vr6a. To leading order, the variance in the image deflection, convergence and shear 
are 

{\5^?)^l^ha {n^) ^ ^l^-ln"- (16) 

l2^c 2 Zjc CL S 

where we must introduce a core radius s to make the variance in the convergence and shear 
finite. These perturbations may then be magnified by the macro model to produce the observed 
perturbations in the image positions and fluxes. If the scale a is determined by the satellite's 
tidal radius, we have a = (bbo)^/'^ where bo is the critical radius of the primary lens (see e.g. 
Metcalf & Madau |2001| ). Thus, for a fixed satellite surface density, the variance of the astrometry 
perturbation in units of the critical radius of the macro lens is roughly 

while the variance in the shear and convergence is roughly 
where InA = ln(-v/66o/s) ~ 10. 
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2.3. Physical Scales and Interpretations 

For this first attempt at modeling substructure, we consider satellites with constant surface 
density S/Sc and critical radius b. We scale the satellite break radius like a tidal radius with 
a = {bboy^'^ for bQ = I'.'O. Near the critical radius of a moderately elliptical isothermal lens the 
surface density is Sc/2, so the projected satellite mass fraction is fsat = 2S/Sc. In the cylinder 
(sphere) defined by the Einstein radius, roughly 50% (10%) of the mass is dark matter (see, e.g., 
Keeton p001b| ). The Einstein ring, where we make the measurement, is typically 1.0-1.5 effective 
radii from the lens center and the dark matter fraction will be significantly higher than these 
average values at the edge of the cylinder where we see the images. Hence, we can interpret fsat as 
the fraction of the dark matter near the Einstein ring in substructure with only modest baryonic 
corrections. 

We use only observational parameters in our models, which means that the physical parameters 
will have small shifts between the lenses we consider because of the changing lens and source 
distances. To provide a sense of the physical scales, consider the parameters for PG1115-I-080 with 
source and lens redshifts of Zg = 1.72 and zi = 0.31 respectively. The inner circular velocity of a 
satellite is Vdrc = 9.7(6/0'.'001)i/2 k^/g^ ^nd its mass is M = 3.4(6/0'.'001)3/2 x IO^/i-^Mq. The 
angles b and o correspond to proper distances of 3(6/0'.'001)/i~^ pc and 100(6/0'.'001)^/^/i~^ pc 
compared to the average distance of 1'.'16 or 3.6/i^^ kpc of the images from the lens center. For the 
other lenses, the distances scale linearly with Dql, the circular velocity scales as {D^s / Dos)^^^"^ 
and the mass scales as Dps/ DqlDls- These changes between lenses are sufficiently small 
compared to our logarithmic uncertainties to ignore. 

We will use satellites with fixed properties in our models, so our estimate of the mass scale 
is a weighted average of the satellite masses. Given our statistical uncertainties models with a 
mass spectrum are unwarranted, and we should be able to estimate the effects of using a mass 
spectrum simply by matching the variance in the shear and astrometry perturbations (Eqn. [l6[) . 
The mass function of the satellite halos is dN/dM oc Af^" with 1.7 ^ a ^ 1.8 (e.g. Moore et 
al. |1999| , Klypin et al. |1999| , Metcalf & Madau |2001| , Springel et al. p001| , Helmi et al. p002| ) which 



we limit to a finite range Miow < M < Mhigh to avoid divergences in the total mass. With a < 2, 
only the upper mass limit is important in estimating the perturbations. Our effective substructure 
mass M is related to the upper mass limit by M = Mhigh(2 — a)/{2> — a) = Mhigh/6 if we match 
the amplitude of the astrometry perturbations and by M = Mhigh((2 — a)/(7/3 — a)Y = Afhigh/20 
if we match the shear perturbations. Given the precision with which we can currently estimate 
the characteristic mass scale, we choose not to include a satellite mass function. Crudely, we can 
estimate that Mhigh ~ 10 - 20M. [| 



^For a = 2 the results become logarithmically sensitive to the lower mass limit. If we match the astrometric 
perturbations, our mass scale corresponds to M = Mhigh/lnC ~ Mhigh/lO where C = Afhigh/Afiow is the ratio of the 
upper and lower mass limits. If we match the shear perturbations, assuming the Coulomb logarithm is held fixed, 
our mass scale corresponds to M = A/high(3/ In C)'^ ~ Mhigh/30. 
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(1) observed lens 
(2) fit SIE + external shear lens model 



(A) Monte Carlo tests 

i 

(A1) add known substructure 

i 

(A2) nonlinearly solve for new images 

i 

(A3) add measurement errors 

i 

(A4) fit SIE + external shear lens model 
(A5) linearized analysis for substructure 
(A6) add to Bayesian likelihood 



(B) final analysis 



V 

(B1) linearized analysis 
for substructure 

t 

(B2) add to Bayesian 
likelihood 



Fig. 1.— 



Outline of analysis procedures. 
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2.4. Analysis Procedures and Monte Carlo Tests 

In Fig. 2 we outline our procedures for analyzing the lens data and for testing our method to 
ensure that it can recover the properties of the substructure accurately by following the treatment 
we use for each lens in the sample. We start (#1) with the available data on the lens. The first 
processing step (#2) is to model the lens with the lensmodel package (Keeton |2001a ) using singular 



isothermal ellipsoids (SIE) in external shear fields for the mass distribution of the primary lens 
galaxy. The SIE model is the only standard lens model which is consistent with general properties 
of the lens sample (see, e.g. Cohn et al. ( p001| ), and references therein). Where needed, we added 
additional SIE lens components so as to reproduce the best standard models for each system. We 
used the observed astrometric uncertainties but broadened the uncertainties in the flux ratios to 
20% to compensate for systematic errors and the contaminating effects of the substructure on 



the flux ratios. Effectively we follow the procedures suggested by Mao & Schneider (1998) for 
modeling lenses in the presence of substructure. 

For analyzing the real data (steps B1-B2 in Fig. |^) we apply our linearized analysis method 
from §2.1 to each lens using the best fit macro model from step ^^2 as the reference model 
(po) supplying the reference image positions and magnifications (x^'^\ M •''''). Given the mass 
scale h and surface number density N of the substructure we determine the angular radius 
Tn = {n/Ny/'^ = {nbaYicrit/^y^'^ inside which we expect to find n perturbing satellites. For each 



realization of the substructure (the 6aj in Eqns. (12) and (13)) we added n = 10 perturbing 
satellites inside radius rio from each image with corrections to avoid over counting in models 
where the rio regions of the individual images overlap. Each satellite perturbed all images, an 
effect which becomes important for mass scales b ^ O'.'Ol. The model is in some senses still a 
"local" approximation because we assume a constant surface density near all images and we do 
not generate a full, global realization of the substructure distribution. The more distant satellites 
produce perturbations which are difficult to distinguish from changes in the macro model. We 
varied only the mass scale b and surface density S/Sc of the substructure using logarithmic 
priors for the two variables {P{b) oc 1/6 and -P(S) tx 1/S). The tidal radius was always set to 
a = (66o)^/^ with bo = I'.'O. For each value of the mass scale and the surface density, we generated 
10^ random realizations of the substructure. All parameters of the macro model are reoptimized 
for every substructure realization by minimizing the fit statistics in Eqns. ( |lO|) and ([Tl| ) combined 
with any ancillary constraints like the position of the primary lens galaxy. The Bayesian likelihood 
distribution (Eqn. 13) is constructed by combining the likelihoods of fitting the data for lens 
j, P{Dj\6aj) for each of the a = 1- ■ ■ 10^ substructure realizations made for each of the set of 
substructure parameters s. 

We tested the algorithm using Monte Carlo models following the steps A 1-6 in Fig. |. The 
objective of the Monte Carlo sequence is to start from the best fit macro model of each lens found 
in step #2 and then by adding substructure and noise generate a synthetic set of lens data which 
should be analogous to the real data. We start by taking the best fit model for the lens (from 
#2) and using its parameters and source properties as the true properties of a new Monte Carlo 
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model. In step Al we randomly place n = 5 perturbing satellites inside the radius of each 
image based on the desired physical properties (mass scale b, surface density S, tidal radius a) of 
the substructure. 



In step A2 we use the lensmodel package (Keeton |2001a ) to find the non-linear solutions for 



the new image positions and fluxes including all the substructure but keeping the macro model 
and source properties fixed. We used n = 5 perturbing satellites per image because of limitations 
on the maximum number of lenses in lensmodel. Tests varying the number of perturbers used 
both to generate and analyze the data suggested that the choices have no effects on the results in 
the sense that any biases are small compared to the statistical uncertainties. A modest fraction 
of realizations for 62045+265 produced extra images. We discarded these realizations. The 
existence of these solutions suggests that the substructure profile shape may be constrained by 
the production of extra images, but an exploration of these additional parameters is beyond the 
scope of our current study. After adding measurement errors to the image positions, lens positions 
and image fluxes in step A3, we have a set of synthetic lens data that should be a realistic Monte 
Carlo model of the real data we started with in step #1. The remainder of the analysis is identical 
to that for the real data. Step A4 matches step ^^2 where we fit the synthetic data using only a 
smooth macro model to provide the reference models and images for performing the substructure 
analysis. Step A5 matches step Bl for the real data, where we apply our linearized substructure 
analysis to the noisy synthetic data and the reference model, and step A6 matches step B2 where 
we combine the results to estimate the Bayesian likelihood distributions for the substructure 
parameters. 

We illustrate the ability of our linearized analysis method to correctly extract the properties of 
substructure in three steps. First, we examined the sensitivity of our surface density estimate fsat 
to measurement errors. Next, we tested our ability to estimate fsat when the satellite masses and 
structures were fixed to their true values. Finally, we tested our ability to estimate simultaneously 
the surface density fsat and the mass scale b. In each case we generate a Monte Carlo data set 
consisting of a perturbed realization for each of the 7 lenses in our sample and then analyze it 
using the same procedures we apply to the real data. 

In our first test we examine whether measurement errors can be mistaken for substructure 
by adding random astrometry and flux errors to the models, fltting new macro models, and 
then analyzing the synthetic data using our method. This is unlikely to be a serious problem 
for the real data because the best fit models typically have ^ ^^doi when we use realistic 
uncertainties for the image fluxes. However, this will determine the range of fsat to which we 
are sensitive, where the relevant scales are 10~^ ^ fsat ^ 10~^ for normal satellite populations 
and 0.02 <^ fsat ^ 0.15 for CDM substructure. The results of two such simulations are shown in 
Fig. ^. The formal, one-sided 90% confldence upper bounds are fsat ^ 0.004, although this is 
very conservative because we only calculated the probability over the range 10^'^ < fsat < 1- The 
peak probability and most of the integrated probability comes from still lower satellite fractions. 
Lenses with highly magnified images are more sensitive to substructure and constrain S/Sc more 
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0.25 




' c sat 



Fig. 2. — Two Monte Carlo simulations of the effects of measurement error. The dashed curves 
show the likelihood distribution for each lens in the first realization and the heavy solid curves show 
the final Bayesian probability distributions for the two realizations of 7 lenses. Our formal upper 
limit is fsat ^ 0.004 in both trials, but the limit would be lower had we extended the calculation 
beyond the range 10"'^ < jsat < 1- 



0.25 



0.2 




/ c sat 

Fig. 3. — Eight Monte Carlo simulations with fsat = 0.05 and b = O'.'OOl. The light dashed 
curves show the likelihood distribution for each lens in the first realization, the heavy solid curves 
show the final Baycsian probability distributions for all eight realizations of 7 lenses each, and the 
heavy dashed curve shows the combined probability of all 8 realizations. This latter case mimics a 
sample of 56 lenses. The points on the heavy curves mark the median probability (triangles) and 
the regions encompassing 68.3% {la, squares), and 95.4% (2a, pentagons) of the likelihood in the 
Bayesian probability distribution. The vertical line marks the true value of fsat = 0.05. 
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strongly, with the upper Umit varying as the inverse square of the maximum image magnification. 
Our lens sample has two "low-magnification" lenses {Mmax ^ 5, B1608+656, B1933+503) and five 
"high-magnification" lenses {M^ax 5, MG0414+0534, B0712+472, PG1115+080, B1422+231, 
and B2045+265). Individual lenses can even show probability peaks at larger surface densities, 
but without the contrast between the peak and the probability at lower fsat needed to produce a 
signal at higher surface density in the full sample. If, however, we underestimate the flux errors, 
we can make a spurious detection of the substructure. When we examine a range of satellite mass 
scales b as well as the surface density, we find that measurement errors produce no preferred mass 
scale for the substructure. 

In the second set of simulations we added perturbing satellites near each image with a surface 
density of fsat = 2S/Ec = 0.05 and a mass profile defined by 6 = O'.'OOl and a = 0'.'032. The radius 
encompassing an average of 5 satellites, = O'.'OSO, is much smaller than the distances between 
the images. We experimented with other models, but the results we present were typical. As 
we show in Fig. |3|, the estimated surface density is consistent with the surface density used to 
generate the data. In the eight simulations shown in Fig. |3|, the surface density corresponding 
to the median probability ranges from fsat = 0.01 to 0.08 with uncertainties of a factor of 2.5 
at la and 3.0 at 90% confidence. The true surface density is within the 68.3% (la) confidence 
region in 4 of the 8 simulations and within the 90% confidence region for 6 of the 8 simulations. 
If we combine all 8 simulations to mimic a sample of 56 lenses, the surface density estimated by 
the median of the Bayesian likelihood distribution is fsat = 0.034 with a 90% confidence range of 
0.023 < fsat ^ 0.048 that marginally excludes the true value. The slightly low value for fsat could 
be due to chance, discarding the cases producing additional images, linearizing the problem or the 
local approximation for the substructure. 

We also examined the likelihood distribution in the two-dimensional space of the surface 
density fsat and the mass scale b, holding the internal structure of the satellites fixed with 
a = (bboy^"^ for bo = I'.'O. Fig. ^ shows that the method can recover the mass scale, but less 
robustly than the surface density of the satellites. In these two-dimensional models we find median 
estimates for the surface density ranging from fsat = 0.014 to 0.074 with a factor of 3.8 uncertainty 
at 90% confidence. The ability to adjust the mass scale significantly increases the uncertainty in 
the surface density. The median estimates for the mass scale range from b = 0'.'00016 to 0'.'0027. 
The uncertainty in the mass scale is usually an order of magnitude at 90% confidence. As we 
discuss in §2.2, we expect to be more sensitive to the surface mass density than to the mass scale. 

—2/3 

For constant astrometry perturbations we expect b oc fsat and for constant shear or convergence 
perturbations we expect b oc f~J- (see eqns. and |l^. Neither slope is clearly reflected in 
the likelihood contours of Fig. ^ suggesting that both types of perturbations contribute. If we 
combine all 8 realizations to mimic a sample of 56 lenses, we recover the input model with modest 
uncertainties. 

In summary, the lenses are sensitive to surface densities of substructure exceeding fsat ^ 0.004 
and samples of 7 lenses can be used to determine the surface density and mass scale with 
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Fig. 4. — Monte Carlo simulations including both the surface density fsat and the mass scale b. 
We show the probability contour enclosing 90% of the total probability for four of the eight models 
shown in Fig. ^. The heavy contour shows the result after combining all eight realizations and the 
solid square marks the true model parameters. 



0.25 



CDM 




Fig. 5. — Results for the observed lens sample with b = (/.'OOl. The heavy solid lines show the 
probability distributions assuming errors in the flux ratios of 5%, 10% and 20%. The points on 
the curves mark the median surface density (triangles) and the regions encompassing 68.3% (la, 
squares), and 95.4% {2a, pentagons) of the probability. The dashed curves show the contributions 
from the individual lenses for the 10% case. The region between the vertical lines is the range 
of substructure mass fractions found in the Klypin et al. ( |199S| ) simulations. Normal satellite 
populations, with 10~^ <i fsat ^ 10^^, correspond to a region off the left edge of the figure. 
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Fig. 6. — Results for the observed lens sample as a function of both the surface density fsat and the 
mass scale b. For the case with 10% flux errors we show the probability contours enclosing 68% (la), 
90%, and 95% {2a) of the total probability using heavy solid curves. For the cases with 5% and 20% 
flux errors we show only the probability contour encompassing 90% of the probability using a light 
solid curve. The region between the vertical lines is the range of substructure mass fractions found 
in the Klypin et al. ( |1999| ) simulations. Normal satellite populations, with 10"^ ^ fsat ^ 10~^, 
correspond to a region off the left edge of the figure. 



-17- 



reasonable accuracy. Our method shows no signs of biases in the recovered parameters given 
the statistical uncertainties expected for 7 lenses. In much larger lens samples, or samples with 
very accurately measured image fluxes, we may underestimate the surface density slightly as a 
consequence of the simplifications made to allow rapid calculation (linearizing the problem and 
the "local" approximation for the substructure) and the elimination of realizations generating 
extra images. 



3. The Properties of Halo Substructure 

Given these limitations we now apply the analysis to the sample of 7 real lenses. We will 
assume an average uncertainty in the flux ratio measurements of 10%, but report results for 
uncertainties of 5% and 20% as well. For the few cases where the flux measurement errors are 
larger we use the measurement errors instead, but in most cases the flux measurement errors are 
dominated by systematic uncertainties rather than measurement errors. Amongst the systematic 
issues are variability and time delays, wavelength dependencies to the flux ratios, and any 
contribution from stellar microlensing. A detailed examination of these problems is beyond the 
scope of our present study. Given the current data for most lenses, the image flux uncertainties are 
certainly lower than 20%, probably lower than 10% and unlikely to be lower than 5%. The errors 
in the image and lens positions are dominated by measurement errors rather than systematic 
errors. 

We first analyzed the data assuming a fixed mass scale of 6 = O'.'OOl and a tidal radius of 
a = 0'.'032. As shown in Fig. ^, the results for the real lens sample have qualitative properties 
that are very similar to the results of the Monte Carlo simulations shown in Fig. The median 
estimate for the surface density depends on the assume level of systematic uncertainties in the 
image flux ratios, with fsat = '^^/'^crit = 0.051, 0.024 and 0.0097 for flux ratio uncertainties of 5%, 
10% and 20% respectively. The 90% confidence ranges for the three cases are 0.027 < fsat < 0.096, 
0.0098 < fsat < 0.058 and 0.0014 < fsat < 0.037 respectively. In all three cases the distributions 
are broadly consistent with the 0.02 < fsat < 0.15 range found in the Klypin et al. ( [L999| ) 
simulations, and well above the 10~^ ^ fsat ^ 10^'^ found in visible satellites (see Mao &; 
Schneider [l998|, Chiba|200lD. 



We also calculated the probabilities as a function of both fsat and the mass scale b as 
shown in Fig. ^. With 10% fiux errors the median estimates for the surface density and mass 
scale are fsat = 0.020 and ^'0013 with 90% confidence regions of 0.0058 < fsat < 0.068 and 
O'.'OOOl < 6 < 0'.'007. For a dn/dM oc l/M^ (Miow < M < Mhigh) satellite mass function this 
implies that the upper mass scale is in a range W^Mq ^ Mhigh ^ W^Mq that is consistent with 
the expectations for satellites. There is a relatively strong covariance between the parameters 
b and fsat, with low surface densities requiring higher mass scales. The slope of the likelihood 
contours is very close to the b oc f~J- slope corresponding to constant shear or convergence 
perturbations (see eqn. rather than the flatter b oc fsat^ slope corresponding to constant 
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astrometry perturbations (see eqn. 17). If we assume 5% flux errors, then the surface density and 
mass scales are restricted to larger values, with 0.013 ^ fsat ^ 0.078 and 0'.'00036 ^ b ^ 0'.'013. 
If we assume 20% flux errors, a broader range is permitted, with 0.0016 ^ fsat ^ 0.051 and 
O'.'OOOOIS ^ b ^ 0'.'0023. These calculations neglect the smoothing effects of the flnite source size 
(A^ ~ 0.01-1.0 milliarcsecond), which will wash out the perturbations from smaller satellites if 
A9 ^ b. With a flnite source size, we would require a larger satellite fraction to produce the same 
perturbations to the images. On a flnal, if qualitative note, the general properties of the likelihood 
distributions for the real data are remarkably similar to those of the Monte Carlo simulations. 



4. Discussion 



CDM simulations generically produce halos in which ~2-15% of the mass is comprised by 
substructure, which is 50-100 times more mass than is observed in the satellites of the Local Group 
(e.g. Moore et al. 1999 , Klypin et al. 1999| ). This substructure problem, and possible conflicts 
between rotation curves and density cusps and the observed and predicted angular momentum 
distributions in spiral galaxies have been interpreted as requiring signiflcant modiflcations to the 
CDM paradigm (e.g. Spergel & Steinhardt [205o| , Bode et al. |00|, Cohn et al |2000D . 



Here we show that the anomalous flux ratios observed in a sample of 7 gravitational lenses 
can be interpreted as requiring a mass fraction of 0.006 < fsat < 0.07 (90% confldence) in satellite 
halos that is remarkably consistent with the CDM predictions. This estimate assumed 10% errors 
(measurement + systematic) in the estimates of image fluxes, but the predicted surface density 
remains consistent with the expectations for CDM over the plausible 5-20% range for these 
uncertainties. The estimates are always well above the 10~^ ^ fsat ^ 10^^ predicted for known 
satellite populations (see Mao & Schneider 1998, Chiba 20011) . This can be consistent with CDM 
and the lower density of Galactic satellites if star formation is suppressed in most such satellites 
as aheady discussed by Klypin et al. (|1999| ) and Bullock et al. (|200q ). For the dn/dM oc M'^-^ 
(Miow < M < Mhigh) mass function expected for satellites (e.g. Moore et al. |1999| , Klypin et 
al. |1999| ) our test provides a rough estimate of the upper mass scale Mhigh — Mq. While 

this is uncomfortably close to the masses capable of disrupting stellar disks and globular clusters 
(e.g. Moore et al. |1999| ), Font et al. ( |2001| ) find that the expected CDM substructure is consistent 
with the survival of thin galactic disks. Thus, our result confirms a surprising if generic prediction 
of CDM models and can be regarded as a major success of the CDM model. By the same token, 
alternatives to CDM which aim to suppress small-scale power (warm dark matter) or to destroy 
small satellites (self-interacting dark matter) are accordingly disfavored. 

We believe that three other explanations, systematic errors in the data, unmodeled, coherent 
structures in the lens and stellar microlensing, are unlikely. While there are systematic errors in the 
lens data, the anomalous fiux ratios which drive the detection of substructure are present at levels 
far above the measurement errors and appear in multiple observations at differing wavelengths 
over periods of years. They may be misinterpreted but not eliminated. They are also unlikely 
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to be due to coherent structures in the lens galaxy. While we analyzed the lenses using singular 
isothermal ellipsoids in an external shear for the macro model, Metcalf & Zhao ( |2002| ) have shown 
that the flux ratios cannot be explained by a broad range of macro models. The typical lens 
galaxy, including all seven discussed here, is an early-type galaxy whose surface brightness profile 
is well modeled by a smooth, elliptical de Vaucouleurs profile (e.g. Lehar et al. 2000| , Kochanek et 



al. 2000| ) with no obvious photometric residuals. Coherent features in the lenses like spiral arms 



would be trivially detected in most cases. Moreover, if we need fsat ~ 0.01 in compact components 
like satellites to perturb the images, we would require a far bigger mass fraction in large scale, 
coherent structures that cannot produce perturbations isolated to a single image. 

The most problematic alternative explanation is stellar microlensing, which is the same 
physical phenomenon but produced by the stellar populations we know to be present in the lens 
galaxy. The basic argument against microlensing is that it has too small a characteristic angular 
scale (^as=microarcseconds) to produce large, long-lived flux ratio anomalies given the sizes of 
typical radio sources. The Compton limit, and direct VLBI observations the lenses, mean that 
typical sources are resolved on scales of 10-1000//as that are large enough to suppress the effects of 
stellar microlensing. The one apparent case of microlensing a radio source, B1600-I-434, is probably 
due to a superluminal sub-component of the radio source where Doppler boosting gives the source 
a smaller effective size and a rapid modulation time scale (see Koopmans & de Bruyn 2000| ) . Even 



in B1600+434, microlensing provides only a ~5% rms variation in the fluxes. Moreover, many of 
the radio lenses also have constant flux ratios on long time scales (years) which are difficult to 
reconcile with producing flux ratio anomalies using the stars. Finally, our method provides an 
estimate for the characteristic mass scale which is grossly inconsistent with stellar microlensing. 
This is reinforced by detailed analyses of B1422-I-231 (Keeton p002 , Bradac et al. 2002| ) which 



find mass scales compatible with CDM substructure but not stellar microlensing. In summary, 
satellites are the most natural explanation, and the required densities are comparable to that 
expected in CDM and higher than that observed in normal satellite populations. The lenses 
cannot address directly whether they are dark or luminous because of the enormous distances. 

Our examination of the problem is a preliminary one, and our estimates can be extended 
and improved if the following points are addressed. First, the entire question of the image fluxes 
and their uncertainties needs to be carefully reconsidered. We used a fixed measurement error of 
10% for the image fluxes, but the estimated surface density and its uncertainties are affected by 
differences between the true errors and the errors used in the analysis. Until now there has been 
little motivation for determining image flux ratios with high precision (say 1% accuracy), but 
improved analyses will need such high precision. Lens monitoring and time delay measurements, 
already important for using the lenses to determine the Hubble constant without the systematic 



problems of the local distance scale (e.g. Schechter 199£), are needed to eliminate the effects 
of source variability on the flux ratios. In optical lenses, observations over a broad range of 
wavelengths are needed to provide accurate corrections for extinction (see Falco et al. 1999| ). 



Second, improved observations of the lenses are needed. The lensed images of the host 
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galaxies of the radio sources, which are relatively easy to observe using deep infrared imaging 



with HST, can be used to constrain the macro model (e.g. Kochanek, Keeton & McLeod 2001). 
Unlike the unresolved images of quasars or the marginally resolved images of radio cores, large 
lensed structures like the host galaxies (A^ ^ O'.'l) constrain the macro model without being 
affected by substructure. Combining the large scale constraints with the compact images allows 
us to probe the substructure while limiting the ability of the macro model to mask its effects. 
Simultaneously, very high resolution, high dynamic range VLBI observations to map thin, extended 
radio structures can be used to extend the search for substructure over larger regions in each lens 



(e.g. Wambsganss & Paczynski |1992| , Metcalf & Madau p001| ). If the VLBI observations can show 



that the anomalous flux ratios are consistent with the geometric structure of the image, then we 
can completely rule out microlensing as an alternative explanation. Finally, careful searches for 
additional, but faint, VLBI images produced by the substructure may be a powerful means to 
constrain the density profiles of the satellites. We have already found that our assumed density 
satellite profile occasionally produces additional, detectable images, suggesting that a shallower 
density profile would be preferred. 

Third, the analysis can be expanded to include complete treatments of the mass spectrum, 
the density profiles of the substructure and the effects of finite sized sources. These additional 
complications were unwarranted in this first calculation because with only 7 lenses all we can 
realistically say we have measured are the average properties of the substructure. Any model 
producing the same average shear and astrometry perturbations should be consistent with 
the data. It is clear from our Monte Carlo simulations, where our model would occasionally 
generate additional images, that the density distribution of the more massive substructures can be 
constrained by limits on the production of extra images. Given our estimated angular scales for 
the substructure perturbations and the dominance of the mass spectrum by the higher mass halos, 
our effects should be little affected by finite sources sizes. If the typical radio source is 1 mas, then 
we are modestly underestimating the surface density. 

Finally, larger samples of lenses can reduce the considerable Poisson uncertainties. At least 
two additional radio quads have been discovered (B0128+437, Phillips et al. |2000| ; and B1555+375 



Marlow et al. |1999|) in the CLASS survey we used as the basis for our analysis, but lack the 



HST imaging data needed to accurately determine the position of the lens galaxy. Two-image 
lenses, while less optimal because of their lower average magnifications, can be included in the 
analysis when additional lensed structures like the images of the quasar host galaxy or VLBI 
subcomponents provide the constraints needed to break the degeneracies between the macro model 
and the substructure we expect for a simple two-image lens. 

Lastly, we note that other probes of substructure may be possible in Local Group galaxies. 
Very recently, Ibata et al. ( |2001a , 2001b) have suggested that the paucity of tidal streamers in the 



Milky Way halo may betray the presence of halo substructure. Johnston et al. ( 20011 ) similarly 



analyze tidal debris from the disrupted Sagittarius dwarf, and find that stars in these tidal tails 
appear to be more scattered than expected for debris orbiting in a smooth halo. Thus there are 
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tantalizing hints of evidence for substructure within our own halo, and further work along this 
avenue may lead to more definite conclusions than is presently possible. Whatever the outcome of 
these local studies, however, only gravitational lenses can detect directly CDM satellites in which 
star formation was suppressed. 
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